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Abstract 

The scaling of the Schrodinger functional coupling is studied numer- 
ically and perturbatively for an SU(3) lattice gauge field coupled to an 
0(a) improved bosonic spinor field. This corresponds to QCD with mi- 
nus two light flavours and is used as a numerically less costly test case 
for real QCD. A suitable algorithm is developed, and the influence of the 
matter fields on the continuum limit and the lattice artefacts are studied 
in detail. 

1 Introduction 

The strong coupling constant a s of QCD is of particular theoretical interest. 
On the one hand, it can be extracted from jet events which are a property 
of the strong interaction at large energies. On the other hand, as has been 
discussed in detail in JjJ , the running coupling can be computed in lattice gauge 
theory. There, the parameters may be fixed in the non perturbative hadronic 
regime, taking as experimental input for example the pion decay constant and 
the masses of the ir, K, D and B. The computation of the running of the 
coupling up to large energies thus provides a quantitative test of the theory, 
which is believed to be fundamental in the hadronic as well as in the high 
energy regime. Furthermore, it is interesting to find out at which energies the 
perturbative behaviour of a given coupling sets in. 

The basic strategy for such a computation has been proposed by Luscher, 
Weisz and Wolff They use a non perturbative definition of the coupling, 
which runs with the spacetime volume. Its evolution is mapped out by a re- 
cursive finite size scaling technique up to large energies where contact with the 
minimal subtraction scheme is made by perturbation theory. The central object 



in this computation is the step scaling function, which can be understood as a 
beta function for finite scale transformations. At each step in the recursive evolu- 
tion of the coupling, the continuum limit is taken. Other key ingredients include 
0(a) improvement and Schrodinger functional boundary conditions || [|, 9. 

This strategy applies to any asymptotically free theory, and it has first been 
tested for the nonlinear 0(3) model in two dimensions pure SU(2) gauge 
theory ||, and pure SU(3) gauge theory R ||, which can also be interpreted 
as the quenched approximation of QCD. In 0], the ALPHA collaboration has 
just published their first quantitative results for the evolution of the coupling 
in QCD with two flavours. 

However, since simulations in full QCD are notorious for high computational 
cost, the data at Nf — 2 do not yet reach very close to the continuum limit in the 
individual steps of the computation. Therefore, we have decided to also study 
the approach to the continuum in a simpler model that is more easily accessible 
to simulation. It differs from the quenched approximation by depending on the 
fermionic determinant. In this computation, the focus is not so much on the 
running coupling as a function of the energy scale, but rather on the details of 
the approach of the continuum limit in the perturbative regime and at slightly 
lower energy. 

The model we investigate in this paper may be viewed as arising from an 
analytic continuation of the flavour number to negative values and in particular 
to Nf — —2. Since a negative power of the fermionic determinant may be 
represented by bosonic spinor fields with the same indices as fermionic fields, 
the name bermions was coined for such theories ]lO| . The main virtue of these 
models is that the interaction term becomes local and thus numerical simulations 
are considerably cheaper than in full QCD. 

In the literature (e.g. O, ^0|), they were mainly considered from an algo- 
rithmic point of view with the idea in mind to extrapolate in Nf from negative 
values to Nf = 2. However, this is problematical, since fermionic zero modes 
may be encountered (for example at small quark masses or in large physical 
volume) , which dominate the dynamics in the theories at negative flavour num- 
bers. Thus, in our work, no extrapolation of results from negative to positive 
values of Nf is attempted or aimed at. 

In |l2| , two of the authors have published results for the step scaling function 
for unimproved Wilson bermions (Nf — —2) in the perturbative regime. Lattice 
artefacts turned out to be very large. As for the quenched approximation and 
for full QCD with two flavours, we now study the 0(a) improved Nf = —2 
theory. The inclusion of the clover term into the bermionic action poses certain 
algorithmic problems that are dealt with in this article. We present a detailed 
study of the performance of the algorithm used for our Monte Carlo simulations. 

An important input for the understanding of our Monte Carlo data comes 
from lattice perturbation theory. The cutoff effects of the step scaling function 
can be computed perturbatively. They are used in the data analysis. The cutoff 
effects have already been estimated in Jl3j to 2-loop order. However, in that 
calculation, the continuum value of the critical mass has been used as a first 
estimate instead of the finite lattice value, which was not yet available (see 
discussion in jl3| ) . The computation of this critical mass is technically more 
involved due to extra tadpole diagrams that emerge from the non vanishing 
background field. Here we present a computation that includes all diagrams 
and completes the study of p3[ . 
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This article is organized in the following way. In the next section, we reflect 
the most important definitions that occurred in our previous articles. In sec- 
tion 3, we discuss the perturbative expansion of the lattice artefacts of the step 
scaling function. After that, the bermion model and the algorithm used in our 
non-perturbative calculations is discussed. In the last section, our numerical 
results are summarized. 



2 Lattice theory 

Since this work extends our earlier work reported in [|[ fj], |l2) , we only briefly 
summarize the necessary notations. For unexplained conventions, we refer in 
particular to {yj, Hi • 

The theory is set up on a four dimensional hyper-cubic Euclidean lattice 
with lattice spacing a and size T x L 3 , T = L, L being an integer multiple of a. 
The gauge field on this lattice is represented by an SU(3) matrix U(x,fi) that 
is defined on every link between nearest neighbour sites x and x + a/2 of the 
lattice (ft denotes the unit vector in the direction fi = 0,1,2,3). Furthermore, 
on the lattice sites reside N{ flavours of mass degenerate fermionic quark fields 
ip(x), which also carry Dirac and colour indices. We do not specify Nf at the 
moment. Later, we want to consider the theory in which Nf is continued to 
negative numbers. This has to be done after the integration over the quark 
fields has been performed. 

The spatial sub-lattices at fixed times Xq are thought to be wrapped on a 
torus. The gauge field thus fulfils periodic boundary conditions in the space 
directions while the quark fields obey periodic boundary conditions in these 
directions up to a phase e l ° pq |. In the time direction, we impose Dirichlct 
boundary conditions. The gauge field at the boundary takes the form 

U(x, fc)Uo=o = exp(aC), 

U(x,k)\ X0=T = exp(aC'). (1) 

The constant diagonal fields C and C can be chosen such that a constant colour 
electric background field is enforced on the system (3) . The boundary conditions 
for the quark fields are discussed in detail in |L5j , The boundary quark fields 
serve as sources for fermionic correlation functions. They are set to zero after 
differentiation. 

The Schrodinger functional is the partition function of the system, 

Z = e- r = J D[U]D[i>]D[^\ e -«>]. (2) 

It involves an integration over the fields with fixed boundary values at xq = 
and xq = T . For the action, we take the sum 

S[U^] = S g [U] + S f [U^,ip] (3) 

of the 0(a) improved plaquette action 

W = 4£>(p)tr(l-tf(p)) W 
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and the fcrmionic action 



S f [U,$,rl>]=Y i $(x)(D + mo)il>(x). (5) 

X 

Here, D is the 0(a) improved Wilson Dirac operator including the Sheikholeslami- 
Wohlert term [ p~7| multiplied with the improvement coefficient c sw and a bound- 
ary improvement term that goes with c t . Details can be found in jl4], [l5|]. As 
discussed for example in || , the leading cutoff effects from the gauge action can 
be cancelled by adjusting the weights w(p) of the plaquettes at the boundary, 
i.e. one sets 

w(p) = ct(go) (6) 

if p is a time-like plaquette attached to a boundary plane. In all other cases, 
w(p) = 1. 

The improvement coefficient c sw has been computed to 1-loop order of per- 
turbation theory with the result |L8|, |l5| 

c S w(ffo) = 1 + 0.26590(7)g 2 , (7) 

independent of N{ to this order. For Nf = and Nf = 2, c sw has also been 
computed non-perturbatively || [l9| . The results of these simulations can be 
represented in the region < go < 1 in good approximation by the rational 
functions 



c S w (50) 



1 - 0.656ffg - 0.152gg - 0.054gg 
= 1 - 0.922g2 : 



. 1 - 0.454 g g - 0.175g 4 + 0.012g 6 + 0.045 go 8 

Csw( 5 o)U f=2 - l- .720g 2 ' (8) 

The boundary improvement coefficients are only known perturbatively. The 
2-loop value for ct depends (in principle) quadratically on Nf and has the form 

ct(ffo) = l + ( -0.08900(5) + 0.0191410(l)iVf)g 2 

+ ( - 0.0294(3) + 0.002(l)iVf + 0.0000(l)iV f 2 )g^. (9) 

c t is known to 1-loop order, 

£ t (go) = l-0.0180(l)g 2 . (10) 

From the Schrodinger functional, a running coupling may be defined by dif- 
ferentiating with respect to the boundary fields. To obtain a complete definition, 
the diagonal matrices C and C and the direction of the differentiation must be 
specified. Here we differentiate along a curve parametrized by the dimension- 
less parameter r\ at the boundary field " A" of reference , which is favoured 
by practical considerations such as mild cutoff effects. With this choice, the 
induced constant colour electric background field can be represented by 

V(x,n) = e aB ^ x \ (11) 

with 

B o = 0, B k = (x C' + (T-x )C)/T. (12) 
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Now, since V = — 91 gf l Z is a renormalized quantity with the perturbative 
expansion T' = g^ 2 T' a + r' x + . . ., a renormalized coupling with the correct 
normalization is defined as 

(13) 

?7=o 

This coupling can be computed efficiently in numerical simulations as the ex- 
pectation value = ^ ff ^ ■ 

For Nf 7^ 0, the coupling depends not only on the scale L but also on the 
mass mi of the quarks, which we define via the PCAC relation [pp| . To this 
end, the fermionic boundary fields of the Schrodinger functional are used to 
transform this operator relation to an identity that holds up to O(o^) between 
improved fermionic correlation functions on the lattice. In section |3J, this will 
be explained in more detail. 

To define the step scaling function a(u), we set u = g 2 (L) and tune m\(L/a) = 
0. Then we change the length scale by a factor 2 and compute the new coupling 
u' — g 2 (2L). The lattice step scaling function £ at resolution L/a is defined as 

E(«,o/L) = ff(2L)L =g2(i) , mi(i/o)=0 . (14) 

These conditions on g 2 and mi fix the bare parameters of the theory. The 
continuum limit a(u) can be found by an extrapolation in a/L. We expect that 
in the 0(a) improved theory T,(u,a/L) converges to er(ti) with a rate roughly 
proportional (i.e. up to logarithms and higher orders) to (a/L) 2 . 




3 Perturbative computation of the cutoff effects 

The size of the cutoff effects in the step scaling function can be estimated in 
perturbation theory. To this end, the relative deviation of the step scaling 
function from its continuum limit is expanded in powers of u, 

E(u,q/L)-g(u) 

o(u, a L) = — 

a(u) 

= [Sw + S n N { ]u+ [S 2 o + S 2 iNi + S22NI] u 2 + 0(u 3 ). (15) 

It turns out to be quite small at 2-loop level. However, it is still necessary 
to extrapolate the Monte Carlo data to the continuum limit by simulating a 
sequence of lattice pairs with decreasing lattice spacing and fixed coupling u. 
One may use the perturbative expansion of 8(u, a/L) to remove the 0(a) cutoff 
effects up to 2-loop order from the non-perturbative values of the step scaling 
function E(m, a/L). 

The 1-loop coefficients 5ij(a/L), first listed in p8|, are shown in table [l]. The 
2-loop coefficients S2j{a/L) have been estimated in However, the parts of 
8 2 j involving contributions from the quarks contain the 1-loop coefficient of the 
critical bare quark mass m c at which the renormalized quark mass vanishes. 
This zero mass condition has to be specified with the cutoff in place. Thus, in 
the expansion 

m c ^m^+mWg 2 + O(g 4 ), (16) 
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we get a/L dependent expansion coefficients mi l \ which only in the limit 
a/L — ► go over to their continuum values, which is rric — at tree level, 
while the 1-loop value TOc can be found in table 1 of jnj. In |l3|, these contin- 
uum values were used to compute the 2-loop coefficients S2j- So, as the authors 
state, the results presented there can only give a first idea of the size of the 



cutoff effects. To obtain the correct values, we need to compute m 
a/L. 



(1) 



at finite 



L/a 5 



10 



Sn 



<20 



$21 



'22 



4 
■5 
6 
7 
8 
9 
10 
11 
12 



-0.01033 
-0.00625 
-0.00394 
-0.00268 
-0.00194 
-0.00148 
-0.00117 
-0.00095 
-0.00079 



0.00002 
-0.00013 
-0.00014 
-0.00014 
-0.00011 
-0.00009 
-0.00007 
-0.00006 
-0.00005 



-0.00159 
-0.00087 
-0.00055 
-0.00038 
-0.00027 
-0.00020 
-0.00015 
-0.00011 
-0.00009 



-0.00069 
-0.00048 
-0.00033 
-0.00021 
-0.00013 
-0.00010 
-0.00007 
-0.00006 
-0.00005 



0.000724 
0.000411 
0.000199 
0.000102 
0.000058 
0.000038 
0.000026 
0.000020 
0.000016 



Table 1: Perturbative results for 5(u,a/L) up to 2-loop order. 



For the quark mass, we adopt the definition of |Tj] based on the PCAC 
relation. We first introduce the bare correlation functions 

/a(*o) = -^E^(^(^)C(y)75^c(z)\, (17) 

x,y,z * ' 

M*o) = 4^^( p ™ 75 ?° ((z) )' (18) 

x,y,z * ' 

where A a and P a denote the axial current and density and £(x) is the functional 
derivative with respect to the boundary quark fields at xq — 0. Now we can 
define the Xq dependent current quark mass 

+ d )f A (x ) + c A ad^d fp{x ) 
m(X0) = 2frM ' (19) 

where do and 8q are the naive forward and backward derivatives on the lattice. 
As an unrenormalized quark mass, we will use m in the middle of the lattice, 
i.e. 

m (-j) for even T/a, 



!(m(I=2)+m(T±2)) foroddr/a .' ( 2 °) 

The O(a) correction of the axial current is proportional to the improvement 
coefficient ca which is 

c A (g ) - -0.00756(1)^ (21) 

to 1-loop order in perturbation theory p"5| |. Here, we are interested in the critical 
bare quark mass m c at which the renormalizcd quark mass is zero. Since ui\ is 
only renormalized multiplicatively, it is sufficient to require ui\ — in order to 
make the renormalized quark mass vanish. 
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The current quark mass may be expanded in powers of (?q, 
m 1 = mf\m ) + m[ 1 \m )gQ + 0(.9q), 



(22) 



where the expansion coefficients depend on the bare quark mass mo. To set up 
perturbation theory, we consider m,Q also as a series, 



m = rn 0) + m^'g^ + 0(g%), 
and expand mi further as 



At) „2 



mi = m 



r(m o) )+ m«(mr) 



d 



(i) 

TUn — — m 



(0) 



(23) 



.g 2 + O(.g 4 ). (24) 



Therefore, the computation of TOc has to be done in two steps. First we 
compute by requiring 



m£ 0) (mf)) = 0, 



then we can determine mi 1 ^ from 



4 X) H 0) ) 



drriQ 



4 0) H 0) ) 



0. 



(25) 



(26) 



The first step is easily done numerically, the results are shown in table [2| The 
second step mainly amounts to expanding /a and fp up to order g$ and requires 
a slightly bigger effort. 



L/a 


(o) 


(i) 


(i) 


4 


-0.0015131 


-0.26667 


0.0027136 


5 


-0.0016969 


-0.26782 


0.0006933 


G 


-0.0006384 


-0.26984 


0.0001990 


7 


-0.0005761 


-0.26995 


0.0000852 


8 


-0.0003209 


-0.27004 


0.0000451 


9 


-0.0002753 


-0.27005 


0.0000026 


10 


-0.0001835 


-0.27006 


0.0000017 


11 


-0.0001561 


-0.27006 


0.0000011 


12 


-0.0001145 


-0.27007 


0.0000008 



Table 2: Perturbative results for the critical quark mass m c up to 1-loop order. 



The expansion of Ja and fp is outlined in Jig ]. After integrating over the 
quark fields and inserting the contractions of the quark and boundary fields, 
one obtains 



fA,p{x ) 



E l(^{P+TP-U(z~ aO,0)S(z,x) 



x,y,z 



TS&yMy-a^O)- 1 } 



(27) 
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where T = 7075 for f a and T = 75 for fp, and the trace is to be taken over the 
Dirac and colour indices only. P + and P_ are the projectors 

P± = i(l±7o), (28) 

and (. . .)q denotes the gauge field average with a probability density propor- 
tional to 

det(L> + to ) exp{S* G [U]}. (29) 

The correlation functions f a and fp in ( p7j ) contain two quantities that have to 
be expanded in powers of go . One is the quark propagator 

S(x, y) = S^(x, y) + S^\x, y)g + S^(x, y)g 2 + O(g 3 ), (30) 

the other one is the gauge field U, which is expanded around the background 
field V 

t/(x,/i) = V(x,y)e^{g aq il {x)} 

= V{x,y){\+ g aq^(x) + g 2 a 2 q^(x) 2 + 0(g%)). (31) 

At 1-loop order, j 'a and fp now get several contributions: 

1. contributions from the first and second order terms of the link variables 
U(z — a0,0) and U(y — aOjO)" 1 , resulting in diagrams la, lb and 2 of 
figure 0, 

2. contributions from contractions between the first order terms of the link 
variables and the first order terms of the quark propagators, resulting in 
diagrams 3a, 3b, 4a and 4b of figure 

3. contributions from the second order terms of the quark propagators, re- 
sulting in diagrams 5a, 5b, 6a and 6b of figure |l|, and 

4. one contribution from the contraction of the first order terms of both 
quark propagators, resulting in diagram 7 of figure ^. 

As already mentioned in section ||, we need a non-vanishing background 
field in order to define the running coupling. Due to the presence of this field, 
we obtain four more contributions. Taking the average (. . .)g leads not only 
to the diagrams in figure [l], but also to contractions of the first order term 
of the link variables and quark propagators with the first order term of the 
total action, including the gauge fixing and ghost terms needed for perturbation 
theory (see || (|). With zero background field, this first order term vanishes, but 
here it has to be taken into account. These contributions result in the diagrams 
of figure ||. Note that only the diagrams containing closed fermion loops depend 
on the number of flavours Nf. So the only Nf dependent contributions to /a 
and fp at 1-loop order come from diagrams 8a and 8b, whereas diagrams 9a 
and 9b are of opposite sign and thus cancel in the sum. Thus becomes Nf 
dependent, 

m W = TO W + m^Nf. (32) 

In contrast to the case of vanishing background field, the propagators are 
not known analytically, so they have to be computed numerically. Here, we 
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Figure 2: Tadpole diagrams contributing to /a(xo) and fp(xo) at one loop 
order of perturbation theory with non-vanishing background held. The dashed 
line represents the ghost propagator. 

have used the method described in [ gjj . Due to these numerical computations, 
computer time is not negligible. For example, on a 200 MHz Pentium PC, the 
computation of /a and fp at L/a — 16 at 1-loop level took us about 16 hours 
of CPU time. As one has to sum over three momentum components and two 
vertex times, the time needed scales asymptotically with (L/a) 5 . 

The 1-loop correlation functions fj± and fp ^ are found by summing all 
the diagrams. We are now able to compute and thus get the deviation 
S(u,a/L) up to 2-loop order. The results are shown in tables [l] and [|. As 
stated before, the 2-loop coefficients 623 are found to be small. 

4 Monte Carlo simulations with Bermions 
4.1 The bermion model 

In the following, wc will especially be interested in the theory in which the 
number of flavours has been continued to Nf = —2 [l(], [l2|]. This has to 
be done after the integration over the fermion fields has been performed. At 
Nf = — 2 the partition function 

Z = [ D[U]e- s * det{D^D) Nf / 2 (33) 
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can be written as 

Z = j D[U]D[<f>+]D[(j)]e- s »- Sh (34) 

with a local bosonic action 

S b [U,cf>]=a i Y,\(D ( t>)(x)\ 2 . (35) 

X 

Note that, in the actual numerical simulation, we have used the hopping pa- 
rameter representation M of the Dirac operator, which is related to D by 

M = 2K(D + m ), k = (8 + 2am y 1 . (36) 

In a similar way, we can express expectation values of fermionic observables 
at Nf = —2 after integration over the quark fields by expectation values in 
the bermion theory. As an a priori guess consistent with 2-loop perturbation 
theory, we chose the improvement coefficient c sw by linearly extrapolating the 
non perturbative results at Nf — 2 and N{ = 0, 

Csw(ffo)| Arf= _ 2 = 2<w(5o)|jv (=0 - c sw(5o)|jv f=2 - (37) 

This choice guarantees that the observables are smooth functions of the bare 
coupling and an extrapolation to the continuum limit is feasible. Furthermore, 
we have also computed the value for c sw for the most critical parameters used in 
this work along the lines of Q and found good agreement with ( |37| ) , see the next 
subsection. For the other improvement coefficients, we take the perturbative 
results with their explicit respective N{ dependence. 

In the bermion model, the occurrence of Dirac operator zero modes is dy- 
namically enhanced. Thus, in situations in which zero modes (or exceptional 
configurations) are to be expected, such as large physical volumes or large values 
for the bare coupling, these may render the functional integral ill-defined or at 
least hamper its Monte Carlo evaluation. However, we will study this theory 
not too far from the perturbative regime so that these problems do not occur. 
Already in the quenched approximation, exceptional configurations occur and 
invalidate measured fermionic correlation functions. To reach larger couplings, 
one could also here consider a twisted mass term pS , 24 , which we shall however 
not pursue in this paper. 



4.2 The size of c sw 

In perturbation theory, c sw is linear in Nf up to two loops, which motivated us to 
first do a linear extrapolation of existing non-perturbative data. To corroborate 
this further, we have computed c sw also non-perturbatively at Nf = —2 for the 
bare coupling (3 — 8.99, which was used in our simulations, compare table [|. 
The computation was done along the lines of ref p2| , to which we refer for 
unexplained notation and an explanation of the method. 

We have computed the current mass aM and the lattice artefact aAM for 
various values of n at three trial values of c sw . As seen before in || and [^9|, we 
note that AM depends only weakly on the mass M. Therefore, we are satisfied 
with a current mass M that roughly vanishes, thereby introducing a negligible 
error. Our results for the lattice artefact are summarized in table II 
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aM 


aAM 


1.171815 


0.0095(1) 


0.00218(19) 


1.271815 


-0.0003(2) 


0.00066(24) 


1.371815 


0.0006(2) 


-0.00127(21) 



Tabic 3: aAM at three trial values of c sw at (3 = 8.99. 



A linear interpolation of these three points to the improvement condition 
aAM — 0.000277 yields c sw = 1.285(7), which is to be compared with the value 
c sw = 1.271815 used in the simulation. The effect of this difference in c sw on the 
step scaling function can be estimated at 1-loop order of perturbation theory. 
For all lattice sizes, the change in £(u, a/L) is smaller than 4 x 10~ 4 , which is 
negligible compared to the statistical errors, compare table ||. 

4.3 Simulation algorithm 

Unimproved Wilson bermions can be simulated with a hybrid overrelaxation 
algorithm, in which the gauge fields arc generated by a combination of heatbath 
and overrelaxation steps and the bermion fields are generated by overrelaxation 
steps only Q . Since the improved bosonic action depends quadratically on an 
individual gauge link U{x 1 /z) through the clover term, we found it practically 
impossible to generalize finite step size algorithms to simulations of improved 
bermions. A local (link by link) hybrid Monte Carlo algorithm would be feasible 
and experience from pure gauge theory shows that it is worthwhile to consider 
it p5| . However, for the same reason as above, a part of the force would have 
to be recomputed at each step on the trajectory. Thus we expect that such an 
algorithm would be very expensive. Also a global hybrid Monte Carlo algorithm 
would be expensive, of course. 

Therefore, we decided to perform global heatbath steps for the bosonic field 
and local overrelaxation steps with respect to the unimproved action for the 
gauge fields. The clover term is then taken into account in an acceptance step 
in which the local action difference with respect to the full improved action is 
used. The acceptance rate in this step turns out to be large enough to pursue 
this algorithm. As the overrelaxation step can be set up symmetrically, the 
combined update fulfils detailed balance. Together with the heatbath step for 
the bosonic field, our algorithm also satisfies ergodicity. 

In principle, boson fields <f>' with the correct distribution can be generated 
by drawing a random field rj from a Gaussian distribution and then applying 
<fi' = A/ _1 7y. However, this procedure is expensive because it requires to run a 
solver with full accuracy for each update. A method published in p6fl, which is 
based on an approximate inversion followed by an additional acceptance step, 
allows to reduce the cost of this step. 

For the update of the gauge fields, we perform a sweep over the lattice and 
update the links sequentially. By an overrelaxation step with respect to the 
unimproved action we propose a new configuration U' which differs from the 
old configuration only for the link variable U(x,fj,) — > U'(x,fi). Thus, for the 
acceptance step, only the part of the action that depends on this link variable 
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is needed. For the gauge part of the action, the difference 

S g [U'} - S g [U] = -|ReTr (([/'(*,//) - U{x,fx))S^(x,fi)) (38) 
can easily be obtained. Here, S^(x,fi) is the sum of the staples at (x,/i) and 

= 

The bermion contribution can also be computed as a local difference. At 
the beginning of a sweep through the lattice, the Dirac operator is applied on 
the whole lattice, and the auxiliary field ip — M[U]cf> is stored. In each local 
step, ip' = M[U']<fi is computed from ip by modifying it at those lattice sites 
that depend on the link variable U(x, fj,) either through the hopping or through 
the clover term. In order to simplify our notation, we split M = Mi + M2 into 
a term M\ which is diagonal in coordinate space and a term M2 which contains 
nearest-neighbour contributions. Then we need to compute 

Aj$#z) = M 1 [U']<f>(z)-M l [U}cj>(z), 

A^(z) = M 2 [U']<t>{z)-M 2 [U]cl>{z). (39) 

Two lattice sites are affected by the hopping term, 

Ai 2 >(z) - -«A M (1 - 7/1 ) (U'(x, n) - U(x, /x)) 4>{x + afi), 
A^x + afi) - -kX* (1 + 7m ) (U'^x, - U\x, »)) cf>(x). (40) 

Here Ao = l,Xk = e l9 . Fourteen lattice sites are affected by the clover term, 
namely x, x + fx and for all directions v ^ [i x ± v and x + /t ± v. At a; we get 
for example 

A^V(x) = ^KC sw (U'{x, fJ ,)-U(x,fj,)) x 

x a flv U(x + a(i, v)U\x + ai>, ^)[/ t (x, v), (41) 

and similar terms at the other points. Since the update is local, one has to be 
careful in parallelizing the algorithm. In a simulation of a lattice on our SIMD 
machine with only two lattice points per node in any direction, neighbouring 
nodes would modify the ip field at a given point simultaneously through the 
clover term. To avoid this conflict, the local lattice size per node in each direction 
has to be larger or equal to three. 

4.4 Performance 

As a measure for the efficiency of our algorithm, we use the machine dependent 
quantity -M cost focusing on the Schrodinger functional coupling g 2 . It is defined 
as 

M cost = (update time in seconds on machine M) 

x (error of l/ff x (4a/T) (4a/ L) 3 . (42) 

In our case, M cos t refers to the CPU time spent on an 8-node machine with 
APE100 architecture. This performance measure allows us to compare for ex- 
ample with performance data obtained in pTj for full QCD. 
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A further important indicator of the efficiency of the algorithm explained 
above is the acceptance rate of the clover term in the gauge field update. This 
acceptance rate turns out to depend only weakly on the parameters in the 
range of couplings considered here. At g 2 = 0.9793, it is about 76%, while at 
g 2 = 1.5145, the acceptance is roughly 70%. 

The cost of our simulations can in principle be optimized by tuning the 
precision of the solver in the boson field update. We could however only obtain 
a total advantage compared to a full precision solver of roughly 10% on the 
small lattices, with a rather flat minimum. This can mainly be attributed to 
the fact that the time of an update step is dominated by the gauge field update. 
Since the optimal size of the residue has to be scaled down when increasing the 



lattice size 26 , this advantage gets even smaller on larger lattices. Hence, we 
expect for our application that tuning the precision parameter on large lattices 
is more expensive than running with an ad hoc guess. 

The cost at g 2 — 0.9793 for various lattice sizes for the improved and the 
unimproved bermion theory is shown in table H. Obviously improvement of 



L/a 


M cost 


M cost 




improved 


unimproved 


4 


0.061(2) 


0.00535(7) 


5 


0.107(3) 


0.00866(13) 


6 


0.212(6) 


0.0155(2) 


8 


0.457(11) 


0.0319(4) 


10 


0.790(17) 




12 


1.30(3) 


0.0788(12) 



Table 4: Costs for improved bermions in comparison with Wilson bermions 
at u — 0.9793. Note that the last two entries for improved bermions are at 
u w 1.11. 



bermions (with the algorithms explained above) leads to a substantial increase 
in computer time, which we estimate to be a factor 12. In section 5[ it will 
become clear however, that it is still profitable. The data of table |i are also 
shown in figure |[ A linear fit in this plot (that excludes the smallest lattices of 
the improved theory) shows that in the improved as well as in the unimproved 
theory, M cost scales as a -2 5 . 

Comparing with data from [£7| , it turns out that improved bermions in our 
implementation are roughly a factor 10 cheaper than simulations in full QCD 
with two flavours. Furthermore, the scaling with a seems to be slightly better 
for the bermions. On the other hand, we have estimated the additional cost of 
unimproved bermions in comparison to pure gauge theory to be only a factor 3. 



5 Technical details and results 

5.1 Parameters 

We have computed the step scaling function S(u, j/) for the two values of the 
coupling u = 0.9793 and u — 1.5145 at lattice sizes L/a = 4, 5, 6, 8. To this end, 
the bare parameters /3 and k have to be tuned such that these couplings are 
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Figure 3: Costs for the measurement of the coupling u = 0.9793 for Wilson 
bermions and improved bermions. 



15 



reached while the current masses m\(L/ a) vanish. The precision required in the 
tuning of n can be estimated in perturbation theory p8| |. In order to estimate 
the effect of a slight mismatch in the tuning of mi, one defines the derivative of 
E with respect to z = L m\ , 



2 (L)=u, mi (L)=z/L 



= Si(o/L) 



(43) 



Under the assumption that it suffices to approximate T,[ by its universal part 
(valid for L/a — > oo), we obtain 



S'i(O) 



N { d 
Att dz 



Cl,l( Z ) 



0.00957 N { . 



(44) 



z=0 



Here, Ci i(z) as defined and computed in [|l6| has been used. This means for 
example that a tuning of the current mass to zero up to 0.001 on an L/a = 8 
lattice leads to an error in the step scaling function smaller than 0.0002 u 2 
and even less on the smaller lattices. Table |5| summarizes the results of the 
tuning procedure. These results allow to neglect the tuning error for the current 
mass while the error of g 2 {L) is propagated into the step scaling function by 
perturbation theory. In those cases where mi = is displayed in the table, this 
is a result of an interpolation in the best tuning runs, which is justified by very 
small values of \ 2 obtained in the interpolation. 



p 


K 


L/a 


9 2 (L) 


mi (L/a) 


10.3488 


0.131024 


4 


0.9793(19) 


0.00000(31) 


10.5617 


0.130797 


5 


0.9795(21) 


0.00055(13) 


10.7302 


0.130686 


6 


0.9793(11) 


0.00000(5) 


11.0026 


0.130489 


8 


0.9793(14) 


0.00000(6) 


8.3378 


0.132959 


4 


1.5145(23) 


0.00000(28) 


8.5453 


0.132637 


5 


1.5145(17) 


0.00000(7) 


8.70830 


0.132433 


6 


1.5145(33) 


0.00000(4) 


8.99 


0.13209 


8 


1.5145(33) 


0.00066(8) 



Table 5: Parameters and results for the coupling and the mass at L. 



5.2 The numerical simulation 

Most of our numerical simulations were performed on APEIOO/Quadrics parallel 
computers with SIMD architecture and single precision arithmetic. We have 
used machines with up to 256 nodes with an approximate peak performance 
of 50 MFlops per node. Roughly half of the statistics for the simulation at 
L/a = 16 and u = 1.5145 has been accumulated on one crate (128 nodes) of 
an APEmille installation in Zeuthen. Since our program was not yet really 
optimized for APEmille, the advantage is only a factor 3. In our simulations, 
we have made much use of trivial (replica) parallelization. 

The coupling and other inexpensive observables have been measured after 
each update, which corresponds to a bosonic heatbath step followed by an over- 
relaxation step for the gauge fields. The fermionic correlation functions to obtain 
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the current mass mi have been measured only rarely, e.g. every 100th update 
sweep, because the mass does not fluctuate much. We have done up to 16 x 31500 
full updates and measurements of the coupling. The statistical errors of the ob- 
servables have been determined by a direct computation of the autocorrelation 
matrix along the lines of appendix A of |^7j . Typical autocorrelation times for 
the coupling range from 3 to 10 (in units of updates). 

5.3 Discussion of results 

For the non-perturbative computation of the step scaling function, we have 
simulated pairs of lattices with size L and 2L at the same bare parameters in 
the bermion theory. The results of these computations arc listed in table ra. For 



L/a 


g J (L) 


9 A $L) 


mi(2L/a) 


4 


0.9793(19) 


1.1090(28) 


-0.00300(10) 


■5 


0.9795(21) 


1.1079(29) 


-0.00086(5) 


G 


0.9793(11) 


1.1053(30) 


-0.00094(4) 


8 


0.9793(14) 


1.1093(40) 


-0.00025(3) 


4 


1.5145(23) 


1.8734(74) 


-0.00266(12) 


5 


1.5145(17) 


1.8648(82) 


-0.00094(7) 


G 


1.5145(33) 


1.8488(86) 


-0.00070(5) 


8 


1.5145(33) 


1.869(14) 


0.00002(5) 



Table 6: Results for the coupling and the mass at 2L at the parameters defined 
by the given value of g 2 (L). 

the propagation of the statistical error and the mismatch of the tuning results 
for the coupling, we use a perturbative ansatz. Then we obtain the lattice step 
scaling function Yl(u,a/L) that is shown in figure We pass to the continuum 
limit by an extrapolation in a/L with the ansatz 

£(u, a/L) = a{u) (1 + p{u){a/L) 2 ). (45) 

As shown in figure ^, this ansatz works perfectly, i.e. within the error bars no 
linear dependence of the step scaling function on a/L can be detected. The 
results for the continuum step scaling function and the corresponding 2- and 
3-loop values are given in table [fl. The difference between 2- and 3-loop pertur- 



u 


a(u) 


^Wh-loop 


c( u )l3-loop 


0.9793 


1.1063(46) 


1.10435 


1.10691 


1.5145 


1.871(17) 


1.85122 


1.87026 



Table 7: Extrapolated simulation results and perturbation theory for the step 
scaling function at Nf = —2. 

bation theory is thus of the same size as the error of the extrapolated simulation 
results. Within the error bars, both values of the step scaling function <j(u) are 
consistent with perturbation theory. 
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Figure 4: Step scaling function for improved bermions for the couplings u = 
0.9793 and u — 1.5145 with fits linear in (a/L) 2 . Shown is also the extrapolated 
continuum value and the 2- and 3-loop results. 
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The typical size of the 0(a) effects can be estimated by considering the step 
scaling function in the unimproved bermion theory as well. This has been done 
in iQ for u = 0.9793. Figure [| shows these results together with the data 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 



Improved and Wilson bermions 

1.16 - 




0.1 0.2 0.3 

a/L 



Figure 5: Results for the step scaling function at u = 0.9793, together with a 
quadratic ht under the constraint of universality. 

after implementing improvement. Obviously, the linear cutoff effects are quite 
seizable for this observable, of the order of a few percent. The data are fitted 
with a combined fit under the constraint that their continuum limit agrees, that 
means assuming universality. This fit is linear plus quadratic in a/L for the 
Wilson bermion data and quadratic in a/L in the improved data. Although 
the additional input from the Wilson data is included, the joint continuum 
limit ^combined (0.9793) = 1.1059(43) agrees almost completely with the value in 
table ^. A linear plus quadratic fit in a/L of the unimproved data alone would 
have given the continuum result <r un i mprovc d(0.9793) = 1.103(12). 

The uncertainties of these extrapolated values show a remarkable success 
of improvement. Although the total computer time for the improved simula- 
tions was only by a factor 1.7 higher than for the Wilson bermions, the error 
after extrapolation is by a factor 2.7 smaller. Since a large portion of the com- 
putational cost for u{u) comes from the largest lattice, this advantage can be 
attributed to the lattice size needed for a reasonable extrapolation. For Wilson 
bermions, this was L/a = 24, whereas simulations are limited to L/a = 16 for 
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the improved case. Of course, our observation is restricted to one value of the 
coupling; at different values, lattice artefacts may behave differently. It should 
also be noted that in the bermion case, adding the clover term leads to a signif- 
icant performance penalty for the numerical simulation. For typical algorithms 
for the simulation of dynamical fermions, like variants of the Hybrid Monte 
Carlo, the inclusion of the clover term implies a much smaller overhead. Hence, 
the advantage of improvement should be even bigger there. 

In addition to our determination of the step scaling function by extrapolating 
L/a) from Monte Carlo simulations, we have also analysed the approach to 
the continuum limit for data which have the 2-loop perturbative lattice artefacts 
cancelled. To this end, we replace the lattice step scaling used for the fit by the 
corrected values 

S(2) («- a / L ) = i w E f!' Q/ w /m 2 ( 4fi ) 
v 1 ' 1 + 5x{a/L)u + 5 2 {a/L)u 2 K ' 

with S from (|l5|) . The results for these fits are shown in figure ||, together with 
the uncorrected fits shown before in figure ||. Again, we leave out the point 
at L/a — 4 for the fits. As can be seen from the plot, this procedure does not 
visibly reduce remaining lattice artefacts. For the coupling u — 1.5145, the slope 
of the fitted line gets slightly smaller, whereas for u — 0.9793 it remains roughly 
equal. However, if we also include the point at L/a — 4, the lattice artefacts 
have even larger 0(a) 2 effects than for the uncorrected data. Nevertheless, 
their continuum limit agrees within the error bars with the one obtained by 
the procedure used before. It is also consistent with perturbation theory in the 
continuum limit. 

Since the mass mi is tuned to zero on the small lattices, we expect that 
L( mi (2L) - mi(L)) is a pure lattice artefact that vanishes in the continuum 
limit with a rate proportional to (a/L) 2 . This expectation is confirmed by 
figure [7] in which this mass difference is shown as a function of (a/L) 2 . While the 
scaling is perfect for the smaller of the two couplings, there are small deviations 
at it = 1.5145, which however can be attributed to statistical fluctuations. 

6 Conclusions 

In this paper, we have investigated the step scaling function in the 0(a) im- 
proved bermion model by means of extrapolating Monte Carlo data at different 
lattice resolutions to the continuum limit. The results obtained were compared 
with renormalized perturbation theory in the continuum limit and with data 
obtained from simulation of unimproved bermions. They also serve as a guide 
in planning analogous simulations with dynamical quarks. 

It is demonstrated that the implementation of improvement successfully re- 
duces lattice artefacts and allows a fit of the step scaling function £(w, a/L) 
linearly in (a/L) 2 (see figure 5). This raises our confidence that the same ex- 
trapolation procedure can be applied for dynamical fermions. We note how- 
ever, that to our disappointment, the computer time required for our algorithm 
turned out to be only about a factor 10 smaller than for dynamical fermions. 
This means that the lattice sizes that can be reached are not much larger than 
for fermions, in contrast to the situation without improvement. 
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Figure 6: Step scaling function for improved bermions for the couplings u = 
0.9793 and u = 1.5145 with fits linear in (a/L) 2 . The rectangles represent the 
data points obtained from our simulations, whereas the circles represent the 
data corrected by perturbation theory. 
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Figure 7: Check of lattice artefacts in the mass at u = 0.9793 and u = 1.5145. 
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